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ABSTRACT 

The intermittent dissipation of interstellar turbnlence is 
an important energy source in the diffuse ISM. Though on 
average smaller than the heating rates due to cosmic rays 
and the photoelectric effect on dust grains, the turbulent 
cascade can channel large amounts of energy into a rela¬ 
tively small fraction of the gas that consequently undergoes 
significant heating and chemical enrichment. In particular, 
this mechanism has been proposed as a solution to the long¬ 
standing problem of the high abundance of CH^ along dif¬ 
fuse molecular sight lines, which steady-state, low tempera¬ 
ture models under-produce by over an order of magnitude. 
While much work has been done on the structure and chem¬ 
istry of these small-scale dissipation zones, comparatively 
little attention has been paid to relating these zones to the 
properties of the large-scale turbulence. In this paper, we at¬ 
tempt to bridge this gap by estimating the temperature and 
CH^ column density along diffuse molecular sight-lines by 
post-processing 3-dimensional MHD turbulence simulations. 
Assuming reasonable values for the cloud density (hu = 30 
cm“®), size {L — 20 pc), and velocity dispersion ((j„ = 2.3 
km s“^), we find that our computed abundances compare 
well with CH"*" column density observations, as well as with 
observations of emission lines from rotationally excited H 2 
molecules. 

1 INTRODUCTION 

The 011*^ ion is commonly detected along sight lines to¬ 
wards bright O and B stars, with column densities > 10^^ 
cm“^ frequently reported in the literature (e.g. Gredel, van 
Dishoeck & Black 1993; Gredel 1997; Crane, Lambert & 
Sheffer 1995; Weselak et al. 2008; Sheffer et al. 2008). This 
prevalence is puzzling, however, because CH^ is destroyed 
very efficiently by both atomic and molecular hydrogen, and 
the only reaction that can form CH'*’ rapidly, 

C+ -b Ha CH+ -b H AE/k = -4640 K, (1) 

is strongly endothermic and can only proceed at temper¬ 
atures of ~ 1000 K or higher. For this reason, models of 
diffuse interstellar clouds with T < 100 K, like those of 
van Dishoeck & Black (1986), fail dramatically to repro¬ 
duce these high CH'*" columns, despite their success with 
other species. 

Most proposed solutions to this problem have invoked 


an additional energy source to overcome this 4640 K acti¬ 
vation barrier. Possibilities include hydrodynamic (Elitzur 
& Watson 1978, 1980) and magnetohydrodynamic (Draine 
& Katz 1986) shock waves, heating in turbulent boundary 
layers at cloud surfaces (Duley et al. 1992), and partic¬ 
ularly dense photon-dominated regions (PDRs) surround¬ 
ing bright stars (Duley et al. 1992; Sternberg & Dalgarno 

1995) ; for an overview of these mechanisms and some of 
the problems they face confronting observations, see Gredel 
(1997). A particularly promising idea, pioneered by Falgar- 
one & Puget (1995), is that the intermittent dissipation of 
turbulence heats small regions within diffuse clouds to the 
> 1000 K temperatures required for (1) to proceed. Drawing 
on laboratory experiments of unmagnetized, incompressible 
turbulent flows, they calculated that if the velocity disper¬ 
sion in cold, mostly atomic clouds at a scale of 1 pc is 3 
km s“^, then a few percent of the cloud could be heated 
to > 1000 K, a mass fraction sufficient to bring the GH'*' 
abundance in line with observed values (Lambert & Danks 
1986). This result was later found to be consistent with mag¬ 
netized, compressible turbulence simulations as well (Pan & 
Padoan 2009). These pockets of warm gas may also explain 
the observed emission from the first few excited rotational 
states of Ha detected in diffuse gas, which is often too large 
to be explained by UV pumping alone (e.g. Falgarone et al. 
2005a; Goldsmith et al. 2010; Ingalls et al. 2011). 

Models that rely on turbulent heating alone can over¬ 
predict the abundance of other species, such as OH, which is 
already well modeled by cold cloud models (Federman et al. 

1996) . However, in addition to the direct heating effect, tur¬ 
bulence can give rise to net drift velocities between the ionic 
and neutral species in plasmas, enhancing the rates of ion- 
neutral reactions like (1) beyond those expected from the 
kinetic temperature alone (e.g. Draine 1980; Flower, Pineau 
des Forets & Hartquist 1985). Federman et al. (1996) ap¬ 
proximated this effect by computing the rate of reaction (1) 
at the effective temperature Tefi given by: 

T,fi = T+^vl, ( 2 ) 

where /r is the reduced mass of (1) and Vd is the magnitude 
of the ion-neutral drift velocity. They proposed that MHD 
waves with amplitudes ~ 3 km s“^ can enhance the pre¬ 
dicted column densities of CH"'" to the observed values even 
in gas that remains T < 100 K. A similar calculation was 
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made in Spaans (1995), who computed the distribution of Vd 
from an analytic intermittency model. More recently, Shef- 
fer et al. (2008) included this effect in their PDR models, 
finding a similar result. The appeal of these models is that 
they have fewer problems over-producing molecules such as 
OH, which is not formed by an ion-neutral reaction. 

The most successful models include both of these effects 
simultaneously. Joulain et al. (1998) and Godard, Falgarone 
& Pineau Des Forets (2009) treat regions of intense dissi¬ 
pation, termed “Turbulent Dissipation Regions” (TDRs) by 
Godard, Falgarone & Pineau Des Forets, as magnetized vor¬ 
tices, taking their (axisymmetric) velocity profiles from that 
of a Burgers vortex, for which the vorticity as a function of 
radius is 

cj(r) = ujo exp - ( — ) . (3) 

\roJ 

Here, cuo and ro are parameters describing the peak vortic¬ 
ity and characteristic fall-off radius in the vortex. Typically, 
cuo « 6 X lO”^*^ s“^ and ro « 40 AU in Godard, Falgarone 
& Pineau Des Forets (2009). These calculations follow the 
subsequent thermal and chemical evolution of parcels of gas 
trapped inside such a vortex, including both turbulent heat¬ 
ing and ion-neutral drift. Godard, Falgarone & Pineau Des 
Forets (2009) then construct models of entire sight lines by 
assuming they intersect some number of these vortex struc¬ 
tures to account for the observed column density of GH"*". 
These models have had a great deal of success reproducing 
the observed CH^ and excited H 2 columns without overpro¬ 
ducing species such as OH. 

The goal of this paper is to provide a complementary 
approach to the above models, which concentrate on individ¬ 
ual dissipation events. We post-process the gas temperature 
T, drift velocity Vd, and CH’'" abundance cell-by-cell through 
an output of a turbulence simulation from Li et al. (2012) 
that has been scaled to typical diffuse cloud conditions. This 
approach loses some of the detail of the above models, but it 
has the advantage of making fewer simplifying assumptions 
about, for example, the nature of the intermittent structures 
or the number of dissipation events along a line of sight. We 
find that CH^ columns in excess of ~ 10^^ cm“^ are read¬ 
ily obtained. We compare our results against a statistically 
homogenous sample of CH’'"-containing sight lines from We- 
selak et al. (2008) and against observations of rotationally 
excited H 2 , finding good agreement with both. 


2 METHODOLOGY 

In this section, we provide an overview of our calculation, 
including our treatment of the heating and cooling rates, the 
drift velocity, and our calculation of the CH^ abundance. 

2.1 Model Description 

The GH^ ion is believed to form in partially molecular en¬ 
vironments. Indeed, in order for reaction (1) to proceed, at 
least some of the hydrogen must be in the form of H 2 , and 
at least some of the carbon must be in G^. Any plausi¬ 
ble formation mechanism is thus not likely to be effective 
in either the outskirts of molecular clouds with visual ex¬ 
tinction Av <0.1 mag, where the hydrogen is almost all 


Table 1. Standard Physical and Chemical Model Parameters 


nn 

30 cm ^ 

L 

20 pc 

ClD 

2.3 km s“ 

Tm 

65 K 

Fso.m 

35 K 

Urms 

5.2 pG 

x(H) 

0.68 

a:(H2) 

0.16 

x(He) 

0.1 

x(e-) 

1.6 X 10“' 

x(C) 

1.6 X 10-' 

x(0) 

3.2 X 10“' 


atomic, or deep in their interiors at Ay greater than a few, 
where almost all the carbon will be G and/or CO. This gas 
is sometimes referred to as the “dark gas” since it is diffi¬ 
cult to observe (Grenier, Casandjian & Terrier 2005; Wolfire, 
Hollenbach & McKee 2010). 

We therefore model interstellar clouds in which the hy¬ 
drogen has begun to turn to H 2 , but the carbon is still 
primarily in the form of C*". Snow & McCall (2006) clas¬ 
sify such clouds as “diffuse molecular clouds.” We treat 
these regions as cubic boxes with length £ 0 , mean hydro¬ 
gen nucleus number density nn, and one-dimensional veloc¬ 
ity dispersion ctid. The mean mass per hydrogen nucleus is 
fiH = 2.34 X 10“^^ g cm“®. For simplicity, we set the rela¬ 
tive abundances (relative to hydrogen nuclei) of molecular 
hydrogen a:(H 2 ) and ionized carbon a;(C''") to be constant 
across the region. For the former, we adopt a:(H 2 ) = 0.16, 
the mean observed molecular fraction from the sample of 
Weselak et al. (2008), which studied the correlation of CH'*’ 
column density with that of atomic and molecular hydro¬ 
gen. For the latter, we take ®(C^) = 1.6 x 10“'* from Sofia 
et al. (2004). Sofia et al. (2011) find using a different mea¬ 
surement technique that the gas-phase carbon abundance is 
lower than that adopted here by a factor of « 0.43. Since it 
is not clear which measurement is more accurate, we choose 
to adopt the higher value. The effects of a lower C abun¬ 
dance in our model are complex. On one hand, it directly 
reduces the CH'*' formation rate (see Section 2.5), since C’*’ 
is one of the reactants in (1). On the other hand, it decreases 
the cooling rate due to C^ (Section 2.3) and increases our 
estimate for the ion-neutral drift velocity (Section 2.4). The 
overall effect of deceasing our assumed carbon abundance 
by a factor of 0.43 is to increase the estimate for the CH'*’ 
abundance by approximately 30%. 

Although chemical and physical models of diffuse gas 
often assume a constant nn, the density distribution in the 
ISM in fact contains a wide range of fluctuations over many 
orders of magnitude due to the compressive effects of su¬ 
personic turbulence. To treat this, we use the results of 
a 512^ driven, turbulence simulation first published in Li 
et al. (2012). The density, magnetic field, and velocity at 
every point in our model are drawn from the corresponding 
cell in a data dump from this simulation, scaled to physical 
units by the process described below. A color plot of the 
column density through the simulation volume is shown in 
Figure 1. An important caveat to our calculation is that this 
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Figure 1. Logarithm of the total column density Nn through the computational domain along the x direction. The simulation data 
has been scaled such that the mean A^h is 1.83 X 10^^ cm“^. The size of the box is indicated on the x and y axes. 


simulation data is isothermal. We then calculate what the 
temperature would be if the intermittency in the isothermal 
case were the same as if the time-dependent heating and 
cooling effects were followed self-consistently. 


2.2 Scaling to Physical Units 

Simulations of magnetized, isothermal turbulent boxes are 
characterized by two dimensionless numbers: the 3D sonic 
Mach number At = ctsd / Cs and the 3D Alfvenic Mach num¬ 
ber J\4a = ctsd/ua. Here, ctsd is the three-dimensional, 
density-weighted, non-thermal velocity dispersion in the 
box, Cs = -^/fer/rh is the isothermal sound speed, where 
T is the temperature and m the mean mass per particle, 
and va = is the Alfven velocity, where Hmia is 

the root-mean-square magnetic field and p the mean den¬ 
sity. In the simulation considered here, the turbulence was 
driven so as to maintain Af « 10, and the initial AIa was 


In the absence of further constraints, we would be free to 
scale p, (Tsd, c^, firms and the size of the box Iq at will as long 
as the dimensionless ratios M and Ma remained invariant 
(see McKee, Li & Klein (2010) for a more rigorous discussion 
of scaling laws for turbulent box simulations). However, to 
be consistent with observations of diffuse molecular gas, we 
impose several additional constraints. First, we require that 
the gas in the box obey a linewidth-size relation (e.g. McKee 
& Ostriker 2007): 

UlD = CTpcfiprf , (4) 

where fipc is the cloud radius in parsecs, and Upc = 0.72 km 
s“^. The ID non-thermal velocity dispersion ctid is related 
to the 3D value by ctsd = \AaiD , and in applying Equation 
(4), which is meant for approximately spherical clouds, to 
our cubic simulation domain, we identify the cloud radius R 
with £o/2. Second, we require that the mean column density 
of hydrogen nuclei be fixed: 


A^h = riuto = Nohs, 


( 5 ) 
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where iVobs ~ 1-83 x 10^^ cm“^ is the mean total column 
density from Weselak et al. (2008). This column corresponds 
to Av ~ 1, consistent with onr reqnirement that the gas be 
partially molecular. Note that Equations (4) and (5) imply 
that we cannot independently choose nn, £o and uid; choos¬ 
ing a density fixes the box size io, which fixes the velocity 
dispersion through the linewidth-size relation. Numerically: 


to 


UlD 


19.8 

2.3 


30 cm 


riH 
30 cm“ 


pc, 


3 \ 0.5 


TlH 


km s 


( 6 ) 


The dynamical timescale in our model can thns be estimated 
as: 


tdyn — 


8.5 X 10“ 


6 / 30 cm 




- 3 \ 0.5 


nn 


yr- 


(7) 


We assume that properties of the fluid flow (the density, 
velocity, and magnetic field) change on this timescale. We 
show that this is large compared to the thermal and chemical 
timescales below. 

These relations also fix (along with the fact the Ma = 
2.2) the rms magnetic field strength in the box: 


firms — 


6-KfJ.HNobaC^^c 

(1 pc) Ml 


5.3 fj,G. 


( 8 ) 


Crutcher et al. (2010) infer that interstellar fields in gas 
with hn < 300 cm“® are uniformly distributed in strength 
between very low values and 10/rG, so this field is quite 
typical. Note that only one of hn, to, usd, and firms may 
be set independently, with the others following from that 
choice. 

The final remaining dimensional parameter describing 
our turbulence simulation is the isothermal sound speed Cs- 
The sound speed of a gas with a:(H 2 ) = 0.16 and a;(He) = 0.1 
is 


However, because the temperature is an output of our model, 
we cannot set it arbitrarily. We thus compute the tempera¬ 
ture using the process described below for a range of boxes, 
each scaled to a different hn, and select the one for which 
A4 computed using the mass-weighted median temperature 
Tso.m (the T for which half of the mass in cloud is hotter) is 
~ 10. We find that this occurs at nu ~ 30 cm“® and adopt 
that as our hducial density. This is the same value adopted 
in the standard model of Joulain et al. (1998). However, in¬ 
dividual sight lines passing through the simulation volume 
can have mean densities ranging from « 5 cm“® to « 180 
cm“®. The corresponding Tsq.m is ~ 35 K, and the mass- 
weighted mean temperature is Tm ~ 65 K. We summarize 
the physical and chemical parameters describing this model 
in Table 1. 

Our simulation data does not include the effects of self¬ 
gravity. As a consistency check, we compute the virial pa¬ 
rameter ftvir to verify that turbulence indeed dominates self¬ 
gravity. Using the above parameters, the total mass M in the 
box is « 8000 Mq. The virial parameter for a spherical cloud 
of mass M and radius R is Qyir = 5a^Y)R/GM (Bertoldi & 
McKee 1992). Using R = ^o/2, we find Ovir ~ 7.4, which is 



T(K) 

Figure 2. Heating (red) and cooling (blue) rates per unit volume 
versus T for nn = 30 cm“® and the standard chemical abun¬ 
dances shown in Table 1. The solid blue curve shows nnAtot, 
while the dashed, dashed-dotted, and dotted curves are tihAhj, 
nHAD+, and rtnAo, respectively. The solid, dashed, and dashed- 
dotted red lines show the mean values of Txurb, hpE, and Tqp, 
respectively. 


large enough that the effects of self-gravity are indeed neg¬ 
ligible. We can also characterize the relative importance of 
gravity and magnetic fields in the box by computing the 
mass-to-flux ratio relative to critical, = M/M^, where 
Mis> ~ ^I2'K'/G and 4? is the magnetic flux threading the 
cloud. For our adopted parameters, ~ 1.3, so the cloud 
is marginally magnetically supercritical. 


2.3 Heating and Cooling 

The temperature in each cell is set by a balance between 
heating and cooling: 

Exurb + TcR -b TpE = nHAtot(T). (10) 

where Exurb, Ecr, and Ere are the heating rates per unit 
volume due to the dissipation of turbulence, cosmic ray 
ionizations, and the photoelectric effect on dust grains, re¬ 
spectively. nnAtot (T) is the total cooling rate per unit vol¬ 
ume, which we assume is dominated by electronic transi¬ 
tions of C'*' and O and by ro-vibrational transitions of the 
H 2 molecule. Note that throughout this paper, we use A to 
represent the cooling rate coefficient (erg cm® s“®) and A 
for the cooling rate per H nucleus (erg s“®). 

Ecr can be expressed as the product of three factors - 
the total cosmic ray ionization rate per H nucleus ((h (in¬ 
cluding both primary and secondary ionizations), the aver¬ 
age energy deposited into the medium per ionization AQ, 
and riH. Both and AQ are rather uncertain and can vary 
considerably over different Galactic environments. Typical 
values of in dense gas are ~ 1 — 5 x 10“®^ s“® (Dalgarno 
2006), but there is evidence from Hs"*" observations that ((h 
is considerably higher in the diffuse gas under consideration 
here (Dalgarno 2006; Indriolo & McGall 2012). Indriolo & 
McGall (2012) find a mean (Ch of 1.8 x 10~®® s“® in their 
sample of diffuse molecular sight lines, and values as large 
as ~ 1 X 10“®® s“® have been reported in the literature 
(Snow & McGall 2006; Shaw et al. 2008). In this paper, we 
adopt the value 1.8 x 10“®® s“®. For AQ, we use 10 eV, as 
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estimated for diffuse molecular gas from Table 6 of Glass- 
gold, Galli & Padovani (2012), although it is important to 
note that this value can vary by several eV depending on 
the precise physical and chemical conditions in the cloud. 
Gombining these factors, the cosmic ray heating rate is: 


PcR = (AQnn 

(11) 

For FpE, we adopt the expression: 

FpE = 1.3 X nueGo ergs cm“® s~^ (12) 

from Wolfire et al. (2003), where Go is the intensity of FUV 
light in units of the Habing (1968) field and e is a heat¬ 
ing efficiency factor given by Equation (20) of Wolfire et al. 
(2003). For nu = 30 cm“^, T — 100 K, an electron fraction 
of 1.6 X 10“^, and a FUV held of Go = 1.1 (Mathis, Mezger 
& Panagia 1983), e evaluates to 1.8 x 10“^, yielding 

PpE = 7.6 X 10“^^ cm“^ (13) 

V30 cm“'^ / 

The hnal heating process we consider is Pxurb- Dimen¬ 
sional arguments (e.g. Landau & Lifshitz 1959) and numer¬ 
ical simulations (e.g. Stone, Ostriker & Gammie 1998; Mac 
Low 1999) both suggest that the kinetic energy in a tur¬ 
bulent cloud l/2pcr|Ei decays in roughly one crossing time 
to/<73D, so that the volume-averaged turbulent heating rate 
is approximately 


Pxurb 


1 P<^3D 

2 £o 


3.5 X 10" 


nu -3 -1 

3 ^ ergs cm s , 


30 cm 


(14) 


where in the last step we have assumed the scaling given by 
Equations (4) and (5). Locally, however, FTurb exhibits large 
huctuations from place to place, a phenomenon known as in- 
termittency. To calculate the spatial dependence of Fxurb, we 
follow the calculation in Pan & Padoan (2009). To summa¬ 
rize their argument, the work done against the viscous forces 
in a Huid with velocity held v is irreversibly converted into 
heat at rate per unit volume given by: 


rTurb(x) = ^diVj + djVi - ^ (V • v) 5ij^ . 


(15) 


This rate depends on the kinematic viscosity of the huid, v. 
However, in our simulations, which were based on the Euler 
equations for a compressible gas, the viscosity was numeri¬ 
cal in origin, and thus does not have its true microphysical 
value. Instead, we treat v as a, proportionality constant that 
takes whatever value is required so that the volume average 
of Equation (15) equals Equation (14). Once this constant 
has been determined, we can compute rTurb(x) for each cell 
in the simulation. 

Note that the resolution of our simulation data Aa; = 
£o/512 « 8 X 10® AU is signihcantly larger than the dissipa¬ 
tion scale provided by the kinematic viscosity of interstellar 
gas of Id ~ 10 AU (Joulain et al. 1998). If the turbulence 
were allowed to cascade down to these small scales, the dis¬ 
tribution of the heating rate would be more intermittent 
than it is in our simulation (Pan & Padoan 2009). How¬ 
ever, viscous dissipation is not the only process that removes 


energy from the turbulent cascade. Ion-neutral friction be¬ 
comes significant at the much larger scale £ad (see Section 
2.4), and is capable of dissipating most (« 70 % for A4a ~ 1) 
of the energy in the cascade at scales « 10 £ad (Li, Myers & 
McKee 2012), which is comparable to the cell size in our nu¬ 
merical data. Below £ad, a turbulent cascade can re-assert 
itself in the neutrals, allowing some of the energy to be dis¬ 
sipated on the ~ 10 AU scale set by viscosity. While our 
procedure likely under-estimates the intermittency for this 
remaining « 30 % of the energy, it is more accurate than the 
opposite assumption that all the energy in cascade makes it 
to ~ 10 AU scales. Note that this procedure for estimating 
the turbulent heating rate is scaled such that it includes all 
the energy removed from the turbulent cascade by dissipa¬ 
tive processes, whether the physical mechanism is molecular 
viscosity or ambipolar diffusion. We do not separately in¬ 
clude an ambipolar diffusion heating term in Equation (10) 
because doing so would amount to double counting. 

For the G"*" and O cooling coefficients, we adopt the 
formulas given by Wolfire et al. (2003): 

Aof (T) = 3.6 X 10~®^exp(—92 K/T) erg cm® s~^ 

X exp (-228 K/T) erg cm® s"\ (16) 


where we have scaled the overall numerical factors to ac¬ 
count for our fractional abundances of carbon and oxygen 
of 1.6 X 10“"^ and 3.2 x 10“'*, rather than the 1.4 x 10“* and 
3.4 X 10“* used in Wolfire et al. (2003). The cooling rates 
per H nucleus are then Aq+{T) = nHAc+(r) and Ao(T) = 
n-HAo(T). These rates are valid for nn < ricrit — 3000 cm“®. 

The final thermal process we consider is the cooling 
rate due to the H 2 molecule. This coolant is particularly im¬ 
portant in the warm (over a few hundred K) gas in which 
GH^ is expected to form. In the low density limit, the 
H 2 level populations are sub-thermal and the cooling rate 
AH 2 (nH —>■ 0)(T) is a sum over the rates due to collisions 
with H, H 2 , He, and e. For these rates, we use the tables in 
Glover & Abel (2008), assuming an ortho:para ratio of 0.7 
(see section 3.3) and the fractional abundances of H, H 2 , He 
and electrons listed in Table 1. These rates are valid at ar¬ 
bitrarily low gas densities, since the cooling rate coefficients 
themselves are independent of the gas density in this limit. 
At densities high enough for local thermodynamic equilib¬ 
rium to be established, the H 2 cooling rate per H nucleus 
Ah 2. LTE(r) becomes independent of the collision partner 
abundances, and we adopt the cooling rate from Goppola 
et al. (2012). We bridge these two limits following Hollen- 
bach & McKee (1979): 


Ah2 (T) 


Ah2, lte(T) 

1 -I- Ah 2 , LTE(r)/AH 2 (riH —>■ 0)(T) 


(17) 


The total cooling rate is then 


Atot (T) = Acf (T) -b Ao (T) -b a:H 2 Ah 2 (F). (18) 

Equation (10) becomes a non-linear equation for T 
in each cell, which we solve numerically using the 
scipy. optimize .brenth routine from the SciPy software li¬ 
brary (Jones et al. 2001-). The magnitudes of these cooling 
processes are summarized as a function of temperature in 
Figure 2 for our standard model parameters. For compari- 



6 


Myers, McKee, & Li 



Figure 3. Blue - the circles show the mass-weighted distribution 
of Vd divided by its mean value (vd) from the A4 = 3, A4a = 0.67, 
^Ad(^o) ~ 1000 AD simulation. The error bars show the 2c7 
temporal variation in distribution over 2 box crossing times, and 
solid line shows the best-fit log-normal. Green - same, but for the 
j\4 = 3, AIa = 0.67 ideal simulation, with computed from 
Equation (22). The agreement between the two curves is quite 
good over more than 3 standard deviations. 



logv^ [cm s M 

Figure 4. Distribution of logu^^ for our model in physical units. 
The blue circles are the simulation data, and the solid line is a 
best-fit normal distribution with a mean of 4.04 and a standard 
deviation of 0.89. 


son, the (constant) values of the various heating rates are 
also displayed for our standard density of uh = 30 cm“^. 

We estimate the typical cooling time in our model as 
the thermal energy density divided by the cooling rate: 

tcooi = 3.3 X loVr, (19) 


2.4 Ion-Neutral Drift 


In the presence of ambipolar diffusion (the net slippage be¬ 
tween the charged and neutral species in a plasma), ion- 
neutral reactions like (1) can proceed at rates faster than 
those expected from the kinetic temperature alone (e.g. 
Draine 1980; Flower, Pineau des Forets & Hartquist 1985). 
The relative importance of the ambipolar and inertial forces 
in a turbulent system can be characterized by the ambipolar 
diffusion Reynolds number i?AD(fo) (Zweibel & Branden¬ 
burg 1997; Li, Myers & McKee 2012): 

( 20 ) 

-^rms 

where pi and are the densities of the ionic and neutral 
components of the fluid and 7 ad is the ion-neutral coupling 
constant given by {av)/{mi + m„). For C*" and H 2 , this 
evaluates to 8.47x10^^ cm® s“^ g“^ (Draine 1980). Applying 
Equations (4) and (5) for our adopted degree of ionization 
and magnetic held strength, i?AD(fo) is 


( 21 ) 


The corresponding length scale ^ad = (-o/R au{D at which 
ambipolar dissipation becomes signihcant is « 640 AU for 
riH = 30 cm“®. Thus, ambipolar drift should not be sig¬ 
nihcant on large scales in our model. However, as with the 
turbulent heating rate, there may be isolated regions in the 
tails of the drift velocity distribution where this effect is 
signihcant. 

To proceed, we need a prescription for computing Vd- 
Unfortunately, two-huid simulations of MHD turbulence are 
prohibitively expensive in the high M, strongly coupled 
regime considered here. To estimate the ehects of Vd on the 
production of CH"*", we instead use our ideal MHD data 
along with an approximate analytic expression for the drift 
velocity in the strongly coupled regime, which we corrobo¬ 
rate with direct numerical simulations of turbulent ambipo¬ 
lar diffusion at lower M . Specihcally, if the system is weakly 
ionized, then the Lorentz force and the ion-neutral drag force 
dominate all the other terms in the ion momentum equa¬ 
tion and the drift is given by Ud = | (V x H) x i3|/47r7ADPiPn 
(e.g. Shu 1992). If the ehects of ambipolar dihusion are weak 
enough that they have only a minor ehect on the geometry 
of the magnetic field, we can estimate the drift by computing 
|(V X H) X H| in the ideal limit: 


Vd 


|(V X H) X 
47r7ADPiPn 


( 22 ) 


This procedure is illustrated in Figure 3. We take two 
simulations of Ad = 3, A4a = 0.67 turbulence from Li et al. 
(2008), one which follows the ion and neutral fluids sepa¬ 
rately, and one which assumes ideal MHD. The AD simula¬ 
tion has i?AD(^o) ~ 1000. We directly compute the time- 
averaged, density-weighted distribution of Vd in the non¬ 
ideal simulation, and compare it to that of Equation (22) 
computed using the ideal data with 7 ad and Xi chosen to 
match the AD simulation. The resulting distributions both 
have an approximately log-normal form: 


where the numerical evaluation is for our fiducial values of 
riH = 30 cm“® and Tsq.m = 35 K. This is two orders of 
magnitude smaller than the dynamical time (Equation (7)), 
which justifies our use of an energy balance equation. 


P{\ogVd)d\ogVd = 

1 

^log v^27r 


exp 


(logUj plogVd) 

2'^log Vd 


(23) 
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and agree with each other to within the error bars, which 
show the magnitude of the temporal fluctuation in the drift 
distribution computed over 2 crossing times. Equation (22) 
does slightly over-predict the simulated value of Vd'- both the 
mean Uiosud standard deviation criog^^ of the Lorentz 
drift approximation are larger than those of the log-normal 
fit to the true drift distribution by 5% and 2%, respectively. 
This is likely because, although i?AD(^o) for the entire box is 
~ 1000 , there are still sub-regions in the box where the cou¬ 
pling less strong. We expect that this approximation should 
improve with increasing i?AD(^o). 

Figure 4 shows the result of applying this procedure to 
our A4 ~ 10, Ma ~ 2.2 Ideal MHD data, scaled to physical 
units as described above. Here, we again And that the dis¬ 
tribution of log(ud) is approximately normal, with best-fit 
parameters uiogv^ ~ 4.04 and uiogv^ ~ 0.89. Although the 
error in the best-fit lognormal parameters found above is 
small, it could potentially have a large impact on the CH"*" 
abundance, since that is primarily determined by the tails 
of the distribution. We can estimate the accuracy of our ap¬ 
proximation as follows. First, we compute the mean CH"*" 
abundance (see Section 2.5 below) in a cloud with constant 
density nn = 30 cm“® and kinetic temperature T = 35 K 
using Equations (26) and (2), under the assumption that the 
distribution of Vd is given by Equation (23) with our best fit 
parameters. We then recompute the CH'*' abundance using 
values of /iiogu^j and criog^^ that are lower by 5% and 2%, 
respectively - the error we found for the A4 = 3, A4a = 0.67 
case above. The result is that the two abundances agree to 
within a factor of 2. As the i?AD(^o) in our target system is 
larger than 1000 (Equation 21), we expect the true error to 
be somewhat less than this. 


2.5 Chemistry 

To assess the viability of turbulent dissipation as an energy 
source for CH"*" production, we perform a simple analytic 
estimate of the CH'*’ abundance in a cell as a function of 
riH and Tefi following Lambert & Banks (1986). CH'*' forms 
through reaction 1 with a rate constant: 

kf = 1.5 X 10“^° X exp(-4640 K/Teff). (24) 

Once it forms, it will be quickly destroyed by reactions with 
H, H 2 , and electrons. The rates for these process are 

CH+-hHC+-hHa, fern = 1.5 X 10"“ cm® s“^ 

CH+ -f Hz ^ cut -h H, fcHa = 1-2 X 10"® cm® s"® 

CH+-h eC-h H, fce = 1.5 X 10~'^ cm® s“®, (25) 

where we have adopted the values used by the Meudon PDR 
code®. 

Because the electron fraction Xe ~ Xi is ^ 10“^, re¬ 
moval of CH"*" by electrons is not crucial and we ignore 
it in our calculations. However, destruction by atomic and 
molecular hydrogen are both important. Balancing the rate 
of formation with the rate of destruction nQ+nns^/ = 


® http://pdr.obspm.fr/PDRcode.html 



Figure 5. Blue solid line - mass-weighted differential distribution 
function of log T. Green solid line - the volume-weighted distribu¬ 
tion of same. Dotted lines - the colors have the same meaning as 
before, but the distribution of logTefj has been plotted instead. 

+ nQjj+tiHzfez 1 derive 

ncH+ + 

« 3.9 X 10-" ( 3 ^) X exp . (26) 

Equation 26 highlights the importance of both the molecular 
fraction and the abundance for CH"*" production; the gas 
must be well-shielded enough from the ambient radiation 
field that some of the hydrogen is molecular form, but not 
so well-shielded that carbon is all molecular. 

We can also estimate a typical CH'*’ formation timescale 
as follows. Suppose we are in a region that initially contains 
no CH^ that is subsequently subject to heating. From reac¬ 
tion 25, the rate of change in is 

—= nc+riHafc/ — nQH+’^Hifcui — (27) 

The solution to this equation is 

ncH+ {t) = ^•CH+,eq X [1 - exp (-at )], (28) 

where eq 4he equilibrium abundance given by Equa¬ 
tion 26 and a = (1 — 22 ;h 2 )fern-I- shz fez- Thus, the time over 
which the CH"® abundance achieves 90% of its equilibrium 
value is 

^ ~ 250 yr. (29) 

This is 2 orders of magnitude shorter than the cooling time 
(Equation 19), and 4 orders shorter than the characteristic 
dynamical time (Equation 7). 

3 RESULTS 

3.1 Temperature and CH^ Abundance 

The result of our temperature calculation is summarized in 
Figures 5 and 6 . Figure 5 shows the mass- and volume- 
weighted differential probability distribution functions of 
logT and logTeff. Figure 6 shows the cumulative distri¬ 
bution functions of the same quantities. Both T and Tea 
can take a large range of values spanning over 3 orders of 
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Figure 6. Same as Figure 5, but showing the cumulative rather 
than the differential distribution functions. Specifically, -P ( > log 
T ) shows the fraction of the mass or volume with values of logT 
greater than the corresponding value on the x-axis. As before, the 
solid lines show the distribution of log T and the dotted lines the 
distribution of logTeff. 



log [cm M 


Figure 8. Blue - average histogram of CH"^ column densities 
over 2500 lines of sight passing through the simulation volume. 
The error bars indicate the 1 — tr range of variation over all 50 
samples of 50 sight lines each. Red - histogram of the CH+ column 
densities from the 50 sight line sample in Weselak et al. (2008). 


magnitude. The distributions of both quantities appear to 
be bimodal, showing peaks at several tens and and several 
hundreds of K. The large majority of the mass lies at low 
T, but a small fraction is found in high temperature “pock¬ 
ets” - precisely the arrangement proposed by Falgarone & 
Puget (1995). An interesting feature of these distributions is 
the importance of density variations in setting the tempera¬ 
ture. The difference between the mass- and volume weighted 
distributions shows that high values of T and Tefr tend to 
be found in relatively low-density gas. This is easily under¬ 
stood from Equations (10) and (22). From Equation (10), 
the heating processes are all proportional to nu. Further¬ 
more, at low T, the dominant cooling processes are and 
O lines, for which the cooling power goes like n^. Accord¬ 
ingly, cooling can balance heating at relatively low temper¬ 
atures unless riH is small. Similarly, the drift velocity Vd is 
inversely proportional to for a fixed ionization fraction 
(Equation 22), meaning that it is likely to be small except 
in low density regions. In our calculations, we find that the 
volume-weighted mean temperature TV « 290 K exceeds the 
density-weighted value Tm ~ 67 K by more than a factor of 
4. Similarly, the volume-weighted median Tso.v ~ 163 K ex¬ 
ceeds the mass-weighted median Tsq.m « 35 K by a similar 
value. This difference is particularly dramatic for the high- 
temperature tails of the distribution: while approximately 
8% of the volume in our box has a Teg greater than 1000 K, 
only 0.3% of the mass does. 

We thus expect CH'^ to mostly be found in low-density 
regions. We show in Figure 7 the regions of T — Vd, T — nn, 
and Vd — nn phase space in which most of the CH^ mass is 
contained. For instance, in the middle panel, we have created 
256^ logarithmically-spaced bins spanning the full range of 
temperatures and densities found in the simulation output. 
The color scale shows the fraction of the total mass in each of 
these T — nn bins, while the black contours show the regions 
of phase space that contain the top 99%, 90%, 50%, 10%, 
and 1% of the the CH^. The left and right panels repeat this 
procedure for T — Vd and Vd — nn, respectively. These plots 


confirm our expectation that regions with T > 1000 K and 
Vd ^ <^30 ~ 4 km s“^ tend to be found at low density, with 
typical values of nn ~ a few H nuclei per cm^. They also 
show that, while these regions are rare, they contain almost 
all of the CH+ molecules. 

The middle panel of Figure 7 illustrates the importance 
of the intermittency in Fxurb in setting the gas temperature. 
The lower curve of this diagram traces out a line that cor¬ 
responds to what the temperature as a function of density 
would be if we considered only Fps, For, and the various 
cooling processes described in section 2.3. Most of the gas in 
the simulation receives little heating from turbulence dissi¬ 
pation and thus lies at or near this line. A small fraction of 
the gas, however, is heated to significantly higher tempera¬ 
tures by the effects of turbulence, and this gas comprises the 
bulk of the material that is heated above ~ 1000 K. We find 
that without FTurb, the fraction of mass in the simulation 
heated to T > 1000 K would drop by approximately a factor 
of 20. 

Sheffer et al. (2008), following Ritchey et al. (2006), 
used the ratio of AfcH+ to A^ch along with Gq and 2 :(H 2 ) 
as an empirical probe of the gas density along a number 
of diffuse sight lines. The resulting density estimates were 
generally quite low, with typical nn ~ 3 cm“®, and in some 
cases much lower than estimates for nn inferred from C*" ex¬ 
citation for the same sight lines (Sonnentrucker et al. 2002, 
2003). Sheffer et al. (2008) interpret this as saying that a 
significant portion of the extinction and atomic hydrogen 
are associated with purely atomic regions that contain no 
CH^, and that the corresponding increase in Go increases 
the estimate for nn- Our result suggests an alternative ex¬ 
planation: that the CH"*" really is predominately found at 
low riH, while the observations are probing something 
closer to the mean density. 

3.2 Sight-line Analysis 

Weselak et al. (2008) presented a sample of 53 CH"*"- 
containing sight lines, 50 of which also had measurements of 
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Mass Fraction [bin ^ ] 

Figure 7. Left - the color scale shows the fraction of the total mass in each logarithmically-spaced T — v^ bin. The black contours show 
the region of phase space in which 99, 90, 50, 10, and 1 percent of the CH+ is found. Middle - same, but for T — riH- Right - same, but 
for Vd-nn. 



Figure 9. CH"^ column density versus total column density A^h 
integrated along the same eight sight lines shown in Table 2. The 
quantities (d) and N-nid] are the CH+ and total column 

densities integrated up though path length d along each ray. The 
x-axis has been normalized by the total column density AIh(^o) 
to fit all the rays on the same plot. 

A^h and A^H 2 ■ For the purpose of comparing with observa¬ 
tions, we will adopt these 50 sight lines as our observational 
sample. We limit ourselves to this sample so that the CH"*" 
column densities would all be calculated in a consistent way, 
but the distribution of columns in Weselak et al. (2008) is 
quite similar to that of Sheffer et al. (2008), which included 
both original data and results from many previous studies 
(Federman et al. 1997; Gredel 1997; Knauth et al. 2001; Pan 
et al. 2004). 

Our model, by construction, has the same mean values 
of Nu and iVna as the Weselak et al. (2008) sample. The 
observed mean and median CH^ column densities for these 
sight lines are 1.2 x 10^^ cm“^ and 1.1 x 10^^ cm~^, respec¬ 
tively. The corresponding values in our model are 1.3 x 10^® 
cm“^ and 8.1 x 10^^ cm“^, i.e. our mean is higher by 8% and 


our median is lower by 26%. To further compare against the 
observational sample, we generate 2500 synthetic observa¬ 
tions by casting 50 groups of 50 rays orthogonally through 
our computational domain as follows. For each ray, we ran¬ 
domly select one of the x, y, and z directions, and then we 
randomly select a coordinate describing ray’s position in the 
corresponding normal plane. For example, for a ray aligned 
with the 2 -axis, we would randomly select a point in the xy 
plane from a uniform distribution to describe the position 
of the ray in the box. For each group of rays, we construct a 
histogram of A^ch+ using 20 bins spaced evenly of the range 
log Vqjj+ = 10 to log VcH+ ~ 15 and compare it against the 
corresponding histogram for the observational sample. The 
result is shown in Figure 8. The error bars indicate the 1 — cr 
variation in the number of sight lines per bin. The number 
of groups was set at 50 because that number was sufficient 
for the error bars to be converged: the maximum change in 
the 1 — a error over all the bins was 15%, and the mean 
change was only 1%. We find that, while our model agrees 
well with the mean and median of the observational sample, 
it does tend to over-produce both very large and very small 
values. However, the sight lines in Weselak et al. (2008) do 
not constitute a truly random sampling in that they were 
chosen because CH+ lines were detectable and the H and 
H 2 column densities were measurable. This is likely not the 
case for the very low column lines in our model. 

From our sample of 2500 rays, we choose 8 “typical” 
sight lines for more detailed investigation as follows. First, 
we randomly draw a ray from the sample. Next, we keep 
it only if both Vch+ and Nn are within 50% of the overall 
sample means. Otherwise, we throw it away and draw again. 
We stop when 8 rays have been selected. The properties of 
these sight lines are described in Table 2. Here, the notation 
xgg means the average value of x in the cells that are in 
the top 99% of the CH'^ number density distribution. Thus, 
h-H,99, 7 m,99, and Vd.m are the mean density, mass-weighted 
temperature, and drift velocity in the regions which contain 
99% of the CH^. We find that, consistent with Figure 7, 
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almost all of the CH"*" is in these sight lines is found in low- 
density, high-temperature, high-drift pockets of gas, with 
typical values of hn.gg « 1 — 2 cm“®, Tivi^gg « 700 K, and 
Vd ,99 ~ 2 — 3 km s“^. These temperatures and velocities 
are quite similar to those obtained in the TDR models of 
Godard, Falgarone & Pineau Des Forets (2009) using very 
different techniques. 

The “temperature” fiuv^/Sk associated with a 4 km s“^ 
drift velocity is « 900 K, comparable to Tm, 99 . Both effects 
thus appear to be important for building up CH'*' columns in 
excess of 10^® cm“^. To gauge the relative importance of the 
two effects, we repeat our calculation with Vd computed as 
above, but with T Hxed at 35 K. The result is that the mean 
CH+ column drops from « 1.3 x 10^^ cm“^ to « 8.0 x 10^^ 
cm“^. If we repeat the same experiment with T computed 
as above, but ignoring the effects of Vd, the CH'^ column 
drops dramatically, down to ~ 5.0 x 10^^ cm“^. Thus, the 
CH+ chemistry in our model appears to be mainly driven 
by ion-neutral drift, with the kinetic temperature making 
a secondary, but not negligible, contribution to the total 
column. 

Finally, to gauge the spatial extent of these regions, we 
show in Figure 9 the result of integrating A^chs- s-nd A^h 
along each of the 8 rays in Table 2. All of the rays show the 
same general behavior: there are large regions that make 
basically no contribution to CH"*" column, punctuated by 
a few thin zones where the CH”*" abundance is substantial. 
The typical sight line intersects approximately 2-4 of these 
regions. This analysis further conhrms the view of Falgar¬ 
one & Puget (1995) - that the cold ISM contains isolated 
patches of hot, chemically active gas, and that these regions 
are crucial to understanding diffuse cloud chemistry. 

The gas densities we infer for the CH+-containing re¬ 
gions in our model are lower than those typically associated 
with gas containing substantial abundances of Hg. As dis¬ 
cussed in Section 2.5, a substantial Hg fraction is crucial 
for CH"*" formation, and our model suggests that low den¬ 
sities are crucial, as well. However, it is not only the local 
gas density that is important for setting the Hg fraction; 
the total shielding from the ambient radiation field mat¬ 
ters as well. In PDR models that assume constant density, 
this distinction is not made, but in a supersonically tur¬ 
bulent medium, it is entirely possible for a region with a 
low local gas density to nonetheless be well-shielded from 
FUV radiation. Indeed, Table 2 and Figure 9 show that, in 
our model, the CH"'’-forming regions tend to be randomly 
distributed along sight lines, so that many of them would 
be well-shielded enough {Ay greater than a few tenths) for 
molecules to form. That said, our adoption of a constant Hg 
fraction is clearly an idealization. A better approach would 
be to self-consistently simulate the formation and destruc¬ 
tion of Hg in the numerical simulation assuming some il¬ 
lumination, so that the molecular fraction (as well as the 
ortho-to-para ratio) could be tracked. 

3.3 Hg Emission 

T > 1000 K, signihcant numbers of Hg molecules can 
be excited to J ^ 2 rotational levels, producing observ¬ 
able emission in the J = 2—>-0, J = 3—>1, and 
J = 4 —> 2 lines. This emission is known to be correlated 
with CH"*" column density (Frisch & Jura 1980; Lambert & 



Ej/k [K] 

Figure 10. Excitation diagram for the first 5 rotational levels 
of Hg from our model compared to observed values. The y-axis 
shows the mean column densities in each rotational level divided 
by the degeneracy factors 3 j. The a:-axis is the energy of each level 
expressed as a temperature. Red circles - Lacour et al. (2005). 
Green circles - Gry et al. (2002). Blue circles - Ingalls et al. (2011). 
Black circles - our model. We have scaled the column densities up 
by a factor of 3 to match the mean Hg column density of the 
observations. The data have been shifted slightly along the r-axis 
for clarity. 


Banks 1986; Jensen et al. 2010), and has been interpreted 
as observational evidence for the intermittent dissipation of 
MHD turbulence (Falgarone et al. 2005b). Because diffuse 
clouds are typically below the critical densities of the J = 3 
and higher lines, the level populations are in general non- 
thermal, and the level populations depend on both T and 
the gas density. Observations of Hg excitation can thus help 
to constrain the temperature and density structure in mod¬ 
els of CH^ production in diffuse clouds. To compare our 
results against these observations, we calculate the H 2 ro¬ 
tational level populations by balancing collisional excitation 
with collisional de-excitation and spontaneous emission for 
the first 20 rotational levels of H 2 . We solve the resulting 
eigenvalue problem using SciPy’s scipy. linalg. eig rou¬ 
tine. The full Cython^ code used for this calculation, as well 
as the rest of the code used in this paper, is available at 
https://bitbucket.org/atmyers/chplus. 

Observations of AfH 2 (-^ = 0) and Ah 2 (J = 1) indicate 
that ortho- (odd J) and para- (even J) hydrogen are not 
generally in found in the equilibrium 3:1 ratio in the ISM. 
For example, in the three sight lines presented in Gry et al. 
(2002), 7 = = 1 )/A'h 2 (J = 0) is 0.7, 0.6, and 0.4. 

Likewise, in the four lines from Lacour et al. (2005), 7 = 0.3, 
1.6, 0.6, and 0.9. Nehme et al. (2008) found 7 = 0.7 towards 
HD 102065, and Ingalls et al. (2011) found that this assump¬ 
tion fit their data for 6 nearby sight lines as well. In this 
paper, we will therefore fix the ortho:para ratio at 0.7 and 
treat ortho-Hg and para-H 2 as separate species. For the en¬ 
ergy levels Ej, we treat both forms of H 2 as quantum rotors 
with rotational temperature Tr = 85.3 K. We used the Ein¬ 
stein A coefficients from Wolniewicz, Simbotin & Dalgarno 
(1998) and the collisional excitation rates from Le Bourlot, 


^ http://cython.org/ 
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iVQH + 

(10^® cm“^) 

A^h 

(1021 cm-2) 

nil 

(cm“^) 

(km s“^) 

7m 

(K) 

UH,99 

(cm“^) 

7m. 99 
(K) 

(km s“^) 

1.1 

2.2 

36.0 

1.8 

60.0 

1.5 

713.0 

2.1 

0.8 

2.2 

36.8 

1.2 

59.5 

1.6 

710.9 

3.0 

1.1 

1.2 

19.8 

1.9 

92.5 

1.7 

708.5 

2.3 

0.9 

1.3 

21.4 

2.2 

85.4 

1.4 

713.0 

2.2 

0.6 

1.9 

31.2 

2.2 

69.3 

1.9 

682.3 

1.7 

1.4 

2.0 

32.4 

0.7 

61.0 

1.3 

715.2 

2.5 

1.6 

1.9 

31.7 

1.8 

66.9 

1.9 

676.6 

2.2 

1.2 

2.7 

43.4 

0.8 

51.0 

1.3 

715.2 

1.9 


Table 2. Data from 8 randomly selected rays cast through the problem domain. Column 1 - the CH+ column density. Column 2 - the 
total column density. Column 3 - the mean number density. Column 4 - The ID rms velocity dispersion. Column 5 - The mass-weighted 
mean temperature. Column 6 - The mean number density in the top 99% of cells by CH+ number density. Column 7 - The mass-weighted 
mean temperature in that same subset of cells. Column 8 - The mean drift velocity in the same cells. 


Pineau des Forets & Flower (1999). The degeneracy factors 
are gj = (2J + 1) (para) and gj = 3(2J + 1) (ortho). 

The result of this calculation for J = 0 to 4 is shown 
in Fignre 10, compared against the H 2 excitation observa¬ 
tions from Gry et al. (2002), Lacour et al. (2005) and Ingalls 
et al. (2011). We have scaled all the column densities in our 
model up by a factor of 3 to match the mean H 2 column 
in the above papers. For all of the levels, our model result 
lies within the spread given by the observations. The mean 
J = 1 to J = 0 rotational temperature in our calculation, 
defined by 


iVH2 (J=l)_ -171 K 

NnAJ = 0) “ Tio 


(30) 


is « 67 K, very close to the mean value for the sample of 
38 sight lines in Rachford et al. (2009) based on FUSE mea¬ 
surements. Because there is a wide range of temperatures 
present and because the J — 2 and higher lines are gen¬ 
erally not in thermal equilibrium, however, our results for 
higher lines are not well-described by a single rotation tem¬ 
perature. 


4 CONCLUSIONS 

The intermittent dissipation of turbulence in the ISM has 
been proposed as an explanation for the high (> 10^^ cm“^) 
CH^ column densities commonly observed along diffuse 
molecular sight lines (Falgarone, Pineau des Forets & Roueff 
1995). Turbulence can aid the production of CH^ both by 
heating a small percentage of the gas directly (Lambert & 
Banks 1986; Falgarone, Pineau des Forets & Roueff 1995; 
Pan & Padoan 2009) and by leading to large drift veloci¬ 
ties (Spaans 1995; Federman et al. 1996; Sheffer et al. 2008) 
within localized regions of intense dissipation. Detailed dy¬ 
namical and chemical models of these intense dissipation 
events have had much success in modeling the chemical 
properties of diffuse sight lines, although the rate of strain 
in the dissipation events and their frequency in the ISM 
were assumed, not self-consistently calculated (Joulain et al. 
1998; Godard, Falgarone & Pineau Des Forets 2009). 

We have re-assessed the origin of the GH"*" observed in 
the diffuse ISM by post-processing a direct numerical sim¬ 
ulation of MHD turbulence, thereby self-consistently deter¬ 


mining the properties of the intermittent turbulence that 
produces the GH’^. We adopted the standard linewidth-size 
relation for molecular gas and set the total column density 
of our simulation equal to the mean value for the observed 
sample that we compare with; the resulting magnetic field 
in our simulation ~ 5 gG is typical for interstellar gas with 
hu < 300 cm“^. We inferred that the density in our simu¬ 
lation is riH — 30 cm“®, which implies that the size of the 
simulated region is £0 ~ 20 pc and the velocity dispersion is 
(JiD ~ 2.3 km s“^. Our approach provides the astrophysical 
framework in which the above models fit, and it corroborates 
their results in several ways. Specifically: 

(i) We have solved an energy balance equation cell-by-cell 
for the temperature. While most of the mass in our cloud is 
cold (< 35 K), a small fraction (< 1%) of the mass has been 
heated to temperatures in excess of 1000 K. 

(ii) Similarly, we have computed the drift velocity in our 
simulation using an approximate analytic expression. While 
on average the drift is negligible, in isolated regions it can 
reach values equal to or exceeding the large scale RMS gas 
velocity of ~ 4 km s“^ . 

(iii) Both of these effects combine to easily produce GH^ 
column densities in excess of 10^^ cm“^. We find that overall, 
the drift is more important, in that it alone accounts for 
about 2/3 of the GH'^ in the box, but the contribution from 
the gas temperature is not negligible. 

(iv) Our work highlights the importance of including den¬ 
sity variations in physical and chemical models of the ISM. 
90% of the GH'^ is found in cells with densities of « 4 cm“^ 
or less - significantly lower than the mean value. These cells 
make up « 5% of the volume and « 0.2% of the mass in 
the simulation and have typical temperatures and drift ve¬ 
locities of « 700 — 800 K and « 3 — 4 km s“^. These values 
are quite similar to the TDR models of Godard, Falgarone 
& Pineau Des Forets (2009), despite the difference of our 
approaches. 

(v) We have estimated the GH^ column density through 
our model and compared the resulting distribution to the 
sample of sight lines presented in Weselak et al. (2008). Our 
mean GH^ column of 1.3 x 10^^ cm“^ agrees very well with 
the Weselak et al. sample, although the median is lower by 
« 26%. 

(vi) Finally, we computed the expected H 2 rotational line 
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emission from these hot regions, and found that it is consis¬ 
tent with observations of diffuse molecular sight lines. 

Finally, there are several caveats to our work that bear 
mentioning. The first, and probably most significant, is 
that our temperature and chemistry calculations were done 
purely in post-processing, taking the magnetic field, veloc¬ 
ity, and density data from an isothermal MHD turbulent 
box. We cannot address whether or not our results would 
have been different had the temperature been computed 
self-consistently and allowed to affect the subsequent flow. 
We have also not self-consistently computed the H 2 fraction 
or the ortho-to-para ratio, instead using an average value 
inferred from observation. However, as chemical models of 
diffuse molecular gas frequently assume constant density, we 
believe that our approach is a valid first step towards a fully 
self-consistent turbulent chemistry calculation. A second, re¬ 
lated caveat is that our results are tied to particular values 
of the sonic and Alfvenic Mach numbers that were used in 
the turbulence simulation. As discussed in Section 2.2, how¬ 
ever, these dimensionless numbers correspond to reasonable 
choices for the cloud density (hu = 30 cm“®), length scale 
(L = 20 pc), and velocity dispersion = 2.3 km s“^). 
Furthermore, our choices for the Mach numbers yield results 
for the mean temperature, CH^ abundance, and H 2 emis¬ 
sion lines that are consistent with observations of diffuse 
molecular regions. A priori, the input parameters needed 
to “predict” these observations could have been inconsis¬ 
tent with the typical properties of diffuse molecular gas, but 
we find that this is not the case. Third, we have used an 
approximate formula to estimate the ion-neutral drift ve¬ 
locity in the strongly-coupled regime (Section 2.4). While 
our formula has been calibrated against low-Mach number 
(A4 ~ 3), non-ideal MHD simulations, we have extrapolated 
these results into the high-Mach number (Af « 10) regime, 
where they have not been directly verified. Doing so would 
require running non-ideal turbulence simulations with am- 
bipolar diffusion in the high-Mach number, strongly coupled 
regime, which are not currently computationally feasible. 
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